Tutorial 4
Jack Halverson
9/27/2021
options(java.parameters = "-Xmx2G")
library(r5r)## Loading required namespace: sf## Loading required namespace: rJava## Please make sure you have already allocated some memory to Java by running:
## options(java.parameters = '-Xmx2G').
## Currently, Java memory is set to -Xmx2Glibrary(osmextract)## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
## Check the package website, https://docs.ropensci.org/osmextract/, for more details.library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.3 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 2.0.1 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(sf)## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1library(ggthemes)
library(ggspatial)
library(tigris)## To enable
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.library(wesanderson)
library(tidytransit)Load Networks
dir.create("networks")## Warning in dir.create("networks"): 'networks' already existsdownload.file("https://app.mecatran.com/urb/ws/feed/c2l0ZT1zbG90cmFuc2l0O2NsaWVudD1zZWxmO2V4cGlyZT07dHlwZT1ndGZzO2tleT0zZTMwMzM1OTRiMTE2NzA0N2IxNjQwNjA0ZjQwMGMzMzdiM2E1MTQ0", file.path("networks","SLOgtfs.zip"), mode = "wb", quiet=TRUE)Street Network Data
SLO_file <- oe_match("San Luis Obispo")## No exact match found for place = San Luis Obispo and provider = geofabrik. Best match is Manitoba.
## Checking the other providers.## An exact string match was found using provider = openstreetmap_fr.Download data to networks folder and read layer of lines representing the street network
SLO_streets <- oe_read(SLO_file$url,
provider = "openstreetmap_fr",
download_directory = "networks",
layer = "lines",
quiet = TRUE) %>%
filter(!is.na(highway))Test Plot
ggplot(SLO_streets) +
geom_sf()
Mask Streets to be only in the City Limits and Project onto Coordinate System
CA5_state_plane <- "+proj=lcc +lat_1=34.03333333333333 +lat_2=35.46666666666667 +lat_0=33.5 +lon_0=-118 +x_0=2000000 +y_0=500000 +ellps=GRS80 +units=m +no_defs"
SLO_city_limits <- places("California") %>%
filter(NAME == "San Luis Obispo") %>%
st_transform(crs = st_crs(SLO_streets))##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|======================= | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|======================================================= | 79%
|
|================================================================== | 95%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%SLO_streets <- SLO_streets[SLO_city_limits,]
ggplot(SLO_streets) +
geom_sf() +
coord_sf(crs = CA5_state_plane)
Load School Locations and rename osm_id column to be “id” for the r5r package to read correctly
SLO_schools <- oe_read(SLO_file$url,
provider = "openstreetmap_fr",
download_directory = "networks",
layer = "points",
quiet = TRUE) %>%
filter(str_detect(other_tags, '"amenity"=>"school"')) %>%
st_filter(SLO_city_limits) %>%
rename(id = osm_id)Create grid of points (hexagons) across entire city to compare travel times to/from
grid <- st_sf(st_make_grid(SLO_city_limits,
square = FALSE,
n = c(100,100),
what = "polygons")) %>%
st_filter(SLO_city_limits)
colnames(grid) <- "geometry"
st_geometry(grid) <- "geometry"
grid <- grid %>%
mutate(id = seq(1, length(grid$geometry), by=1))
grid_points <- st_centroid(grid)## Warning in st_centroid.sf(grid): st_centroid assumes attributes are constant
## over geometries of xggplot(grid) +
geom_sf() +
geom_sf(data = grid_points, alpha = 0.5) +
geom_sf(data = SLO_streets) +
coord_sf(crs = CA5_state_plane) +
theme_map()
Set up r5r core
r5r_core <- setup_r5("networks", verbose = FALSE)## Using cached version from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/r5r/jar/r5-v6.2-all_20210408.jar##
## Using cached network.dat from networks/network.datCalculate travel time matrix
ttm <- travel_time_matrix(r5r_core = r5r_core,
origins = SLO_schools,
destinations = grid_points,
mode = c("WALK", "TRANSIT"),
departure_datetime = as.POSIXct("15-11-2021 14:00:00",
format = "%d-%m-%Y %H:%M:%S"),
max_walk_dist = 1000,
max_trip_duration = 480,
verbose = FALSE)## Warning in assert_points_input(destinations, "destinations"): 'destinations$id'
## forcefully cast to character.Reformat ttm to have each row be an origin and columns for the schools
tt_wide <- ttm %>%
pivot_wider(names_from = fromId,
names_prefix = "from", values_from = travel_time) %>%
rename(id = toId) %>%
merge(grid) %>%
replace(is.na(.), 999) %>%
rowwise() %>%
mutate(from_any = min(c_across(starts_with("from")), na.rm = TRUE))
st_geometry(tt_wide) <- "geometry"Plot travel time to nearest school
ggplot(SLO_streets) +
geom_sf(data = tt_wide,
aes(fill = from_any),
color = NA) +
geom_sf(alpha = 0.5) +
scale_fill_gradient2(low = "green", mid = "yellow", high = "red",
midpoint = 30,
name = "Transit Travel\ntime to the\nnearest school\n(minutes)") +
coord_sf(crs = CA5_state_plane) +
theme_map()
Create Isochrones
iso_pallete <- wes_palette("Zissou1", n = 5)
iso10min <- tt_wide[tt_wide$from_any < 11,] %>%
st_union()
iso20min <- tt_wide[tt_wide$from_any < 21,] %>%
st_union()
iso30min <- tt_wide[tt_wide$from_any < 31,] %>%
st_union()
ggplot(SLO_streets) +
geom_sf(data = iso30min,
aes(fill = "Area within 30 minutes"),
color = NA) +
geom_sf(data = iso20min,
aes(fill = "Area within 20 minutes"),
color = NA) +
geom_sf(data = iso10min,
aes(fill = "Area within 10 minutes"),
color = NA) +
geom_sf(alpha = 0.5) +
scale_fill_manual(values = c(iso_pallete[1],
iso_pallete[3],
iso_pallete[5]),
name = "Transit Travel\ntime to the\nnearest school\n(minutes)") +
coord_sf(crs = CA5_state_plane) +
theme_map()
Load transit stop locations
SLO_transit <- read_gtfs(file.path("networks", "SLOgtfs.zip"))
transit_stops <- st_as_sf(SLO_transit$stops,
coords = c("stop_lon", "stop_lat"),
crs =st_crs(grid))Set up grid of points with stops data
transit_grid <- grid %>%
mutate(num_stops = lengths(st_covers(grid, transit_stops)))
transit_points <- st_centroid(transit_grid)## Warning in st_centroid.sf(transit_grid): st_centroid assumes attributes are
## constant over geometries of xggplot(transit_points) +
geom_sf(aes(color = as.character(num_stops))) +
scale_color_manual(values = c("gray", "cornsilk", "lightgreen"),
name = "Number of\ntransit stops") +
theme_void()
Calculate accessibility - 10 minute walk to transit stop
transit_access <- accessibility(r5r_core,
origins = transit_points,
destinations = transit_points,
mode = "WALK",
opportunities_colname = "num_stops",
decay_function = "step",
cutoffs = 11,
departure_datetime = as.POSIXct("15-11-2021 14:00:00",
format = "%d-%m-%Y %H:%M:%S"),
max_walk_dist = 500,
time_window = 60,
percentiles = 50,
verbose = FALSE) %>%
mutate(id = as.numeric(from_id)) %>%
merge(grid)## Warning in assert_points_input(origins, "origins"): 'origins$id' forcefully cast
## to character.## Warning in assert_points_input(destinations, "destinations"): 'destinations$id'
## forcefully cast to character.st_geometry(transit_access) <- "geometry"
ggplot(transit_access) +
geom_sf(aes(fill = accessibility), color = NA) +
scale_fill_viridis_c(name = "Transit stops\nwithin 10-minutes\nwalk") +
coord_sf(crs = CA5_state_plane) +
theme_void()
Calculate access with a decay function (more continuous) with “accessibility” decreasing by one half every five minutes
transit_access2 <- accessibility(r5r_core,
origins = transit_points,
destinations = transit_points,
mode = "WALK",
opportunities_colname = "num_stops",
decay_function = "exponential",
cutoffs = 5,
departure_datetime = as.POSIXct("15-11-2021 14:00:00",
format = "%d-%m-%Y %H:%M:%S"),
max_walk_dist = 500,
time_window = 60,
percentiles = 50,
verbose = FALSE) %>%
mutate(id = as.numeric(from_id)) %>%
merge(grid)## Warning in assert_points_input(origins, "origins"): 'origins$id' forcefully cast
## to character.## Warning in assert_points_input(destinations, "destinations"): 'destinations$id'
## forcefully cast to character.st_geometry(transit_access2) <- "geometry"
ggplot(transit_access2) +
geom_sf(aes(fill = accessibility), color = NA) +
scale_fill_viridis_c(name = "Accessiblity score") +
coord_sf(crs = CA5_state_plane) +
theme_void()
Turn off r5r core to save memory
stop_r5(r5r_core)## r5r_core has been successfully stopped.rJava::.jgc(R.gc = TRUE)Save grid of accessibility values to use in next assignment
st_write(transit_access2, 'SLO_access.geojson', append=FALSE, quiet=TRUE )## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Error 6: DeleteLayer() not supported by
## this dataset.